R Notebook

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(glue)
library(SingleCellExperiment)

## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics

## Warning: package 'MatrixGenerics' was built under R version 4.3.1

## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: GenomicRanges

## Warning: package 'GenomicRanges' was built under R version 4.3.1

## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: S4Vectors

## Warning: package 'S4Vectors' was built under R version 4.3.1

## 
## Attaching package: 'S4Vectors'
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Loading required package: IRanges

## Warning: package 'IRanges' was built under R version 4.3.1

## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:glue':
## 
##     trim
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## Loading required package: GenomeInfoDb

## Warning: package 'GenomeInfoDb' was built under R version 4.3.1

## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

library(lemur)

## 
## Attaching package: 'lemur'
## 
## The following object is masked from 'package:dplyr':
## 
##     vars
## 
## The following object is masked from 'package:ggplot2':
## 
##     vars

source("util.R")

fit <- qs::qread("../benchmark/output/cable_spatial_plaque_data/fit_small.qs")
nei <- qs::qread("../benchmark/output/cable_spatial_plaque_data/nei_small.qs")

mouse_labels <- structure(paste0("Mouse ", 1:4), names = paste0("mouse_", 1:4))

spatial_plaq_dens_pl <- as_tibble(colData(fit)) %>%
  ggplot(aes(x = x, y = y)) +
    ggrastr::rasterize(geom_point(aes(color = plaque_density), size = 0.05, stroke = 0), dpi = 600) +
    scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1)) +
    facet_wrap(vars(sample), labeller = as_labeller(mouse_labels)) +
    small_axis(label = "spatial coord.", fontsize = font_size_small) +
    labs(color = "Plaque density", subtitle = "Spatial structure") +
    guides(color = guide_colorbar(barheight = unit(1, "mm"))) +
    theme(legend.position = "top", strip.text = element_text(size = font_size_tiny, margin = margin(b = -1, unit = "mm")),
          plot.subtitle = element_text(size = 7))

spatial_plaq_dens_pl

set.seed(1)
umap <- uwot::umap(t(fit$embedding))

as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = sample), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = cell_type_lumped), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = plaque_cluster), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

umap_plt <- as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterize(geom_point(aes(color = plaque_density), size = 0.1, stroke = 0, show.legend = FALSE), dpi = 600) +
    scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1))  +
    small_axis(label = "UMAP", fontsize = font_size_small) +
    labs(subtitle = "Latent structure") +
    theme(plot.subtitle = element_text(size = 7))

umap_plt

Make cell type plot

table(fit$colData$cell_type_lumped, fit$colData$cell_type, useNA = "always")

##                  
##                   ASTROCYTE CHOROID_PLEXUS ENDOTHELIAL_STALK ENDOTHELIAL_TIP
##   ASTROCYTE            8659              0                 0               0
##   MICROGLIA               0              0                 0               0
##   NEUROGENESIS            0              0                 0               0
##   NEURON                  0              0                 0               0
##   OLIGODENDROCYTE         0              0                 0               0
##   Other                   0              2               259             790
##   <NA>                    0              0                 0               0
##                  
##                   MACROPHAGE MICROGLIA MURAL NEUROGENESIS NEURON
##   ASTROCYTE                0         0     0            0      0
##   MICROGLIA                0       816     0            0      0
##   NEUROGENESIS             0         0     0          519      0
##   NEURON                   0         0     0            0  23032
##   OLIGODENDROCYTE          0         0     0            0      0
##   Other                   61         0    86            0      0
##   <NA>                     0         0     0            0      0
##                  
##                   OLIGODENDROCYTE POLYDENDROCYTE  <NA>
##   ASTROCYTE                     0              0     0
##   MICROGLIA                     0              0     0
##   NEUROGENESIS                  0              0     0
##   NEURON                        0              0     0
##   OLIGODENDROCYTE            6308              0     0
##   Other                         0            523     0
##   <NA>                          0              0     0

celltype_broad_plt <- as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterize(geom_point(aes(color = cell_type_lumped), size = 0.1, stroke = 0, show.legend = FALSE), dpi = 600) +
    ggrepel::geom_label_repel(data = . %>% summarize(umap = matrix(colMedians(umap), nrow = 1), .by = cell_type_lumped), 
                              aes(label = cell_type_lumped, color = cell_type_lumped), size = font_size_small / .pt, show.legend = FALSE, max.overlaps = 100) +
    scale_x_continuous(limits = range(umap[,2]) * 1.5) +
    small_axis(label = "UMAP", fontsize = font_size_small) +
    labs(title = "(B) Mouse hypocampus: broad cell type labels") 

cowplot::save_plot("../plots/suppl_alzheimer_cell_types.pdf", celltype_broad_plt, base_height = 8, base_width = 8, units = "cm")

genes_of_interest <- nei$name[1]
mask <- matrix(NA, nrow = length(genes_of_interest), ncol = ncol(fit), 
               dimnames = list(genes_of_interest, colnames(fit)))
mask2 <- matrix(1, nrow = length(genes_of_interest), ncol = ncol(fit), 
               dimnames = list(genes_of_interest, colnames(fit)))
for(id in genes_of_interest){
  mask[id, filter(nei, name == id)$neighborhood[[1]]] <- 1
  mask2[id, filter(nei, name == id)$neighborhood[[1]]] <- NA
}

psce2 <- glmGamPoi::pseudobulk(SingleCellExperiment(list(inside = as.matrix(logcounts(fit)[genes_of_interest,,drop=FALSE] * mask),
                                                        outside = as.matrix(logcounts(fit)[genes_of_interest,,drop=FALSE] * mask2))),
                      group_by = vars(sample, plaque_cluster), n = n(),
                      aggregation_functions = list(.default = \(...) matrixStats::rowMeans2(..., na.rm = TRUE)),
                      col_data = as.data.frame(colData(fit)))

## Aggregating assay 'inside' using 'custom function'.

## Aggregating assay 'outside' using 'custom function'.

comparison_data <- as_tibble(colData(psce2)) %>%
  mutate(expr_inside = as_tibble(t(assay(psce2, "inside"))),
         expr_outside = as_tibble(t(assay(psce2, "outside")))) %>%
  unpack(starts_with("expr"), names_sep = "-") %>%
  pivot_longer(starts_with("expr"), names_sep = "[-_]", names_to = c(".value", "origin", "symbol")) 



expr_comparison_pl <- comparison_data %>%
  mutate(plaque_cluster = ordered(plaque_cluster, levels = levels(fit$colData$plaque_cluster))) %>% 
  ggplot(aes(x = plaque_cluster, y = expr)) +
    geom_point(aes(color = origin), size = 0.3) +
    geom_smooth(aes(color = origin, x = as.integer(plaque_cluster)), span = 1.5, se = FALSE, linewidth = 2) +
    scale_color_manual(values = c("inside" = "black", "outside" = "lightgrey"), labels = c("inside" = "Cells in neighborhood", "outside" = "All other cells")) +
    scale_y_continuous(limits = c(0, 0.3), expand = expansion(add = 0)) +
    guides(x = guide_axis(angle = 45)) +
    theme(axis.title.x = element_blank(), legend.position = "bottom") +
    labs(color = "", y = "Expression",
         subtitle = "\\emph{Jun} expr. vs. plaque density")  +
    theme(plot.subtitle = element_text(size = 7))

expr_comparison_pl

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

is_inside <- tibble(symbol = genes_of_interest, cell_id = list(colnames(fit))) %>%
  left_join(dplyr::select(nei, name, neighborhood), by = c("symbol"= "name")) %>%
  mutate(inside = map2(cell_id, neighborhood, \(ref, nei_cells) ref %in% nei_cells)) %>%
  dplyr::select(-neighborhood) %>%
  unnest(c(inside, cell_id))

de_plot_data <- as_tibble(colData(fit), rownames = "cell_id") %>%
  mutate(umap = umap) %>%
  mutate(de = as_tibble(t(assay(fit[genes_of_interest,], "DE")))) %>%
  unnest(de, names_sep = "-") %>%
  pivot_longer(starts_with("de-"), names_sep = "-", values_to = "de", names_to = c(NA, "symbol")) %>%
  inner_join(is_inside, by = c("symbol", "cell_id")) %>%
  sample_frac(size = 1) 

abs_max <- max(abs(quantile(de_plot_data$de, c(0.95, 0.05))))

jun_de_pl <- de_plot_data %>%
  mutate(inside = ifelse(inside, "in", "out")) %>%
  mutate(inside = factor(inside, levels = c("in", "out")))  %>%
  ggplot(aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterise(geom_point(aes(color = de), size = 0.05, stroke = 0), dpi = 600) +
    # scale_colour_gradient2_rev(limits = c(-1, 1) * abs_max, oob = scales::squish, breaks = c(-1, 0, 1) * signif_to_zero(abs_max, 1), mid = "lightgrey") +
    scale_color_de_gradient(abs_max, mid_width = 0.2) +
    facet_grid(vars(inside), labeller = labeller(inside = as_labeller(c("in" = "Cells in Neighb.", "out" = "All other cells"))),
               switch="y") +
    small_axis("UMAP", fontsize = font_size_small) +
    theme(legend.position = "bottom", legend.margin = margin(t = -2, unit = "mm")) +
    guides(color = guide_colorbar(barheight = unit(1, "mm"), barwidth = unit(15, "mm"), title.vjust = 1)) +
    labs(color = "$\\Delta$")

jun_de_pl

rel_plt <- colData(fit) %>%
  as_tibble(rownames = "cell_id") %>%
  left_join(is_inside) %>%
  dplyr::count(plaque_cluster, symbol, inside) %>%
  mutate(frac = n / sum(n), .by = c(symbol, inside)) %>%
  mutate(plaque_cluster = fct_relabel(plaque_cluster, \(x) str_replace(x, ",", ", "))) %>%
  ggplot(aes(x = plaque_cluster, y = frac)) +
    geom_col(aes(fill = inside), position = position_dodge2()) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
    scale_x_discrete(position = "top") +
    scale_fill_manual(values = c("TRUE" = "black", "FALSE" = "lightgrey")) +
    guides(x = guide_axis(angle = 90)) +
    theme(axis.title.x = element_blank(), legend.position = "bottom") +
    labs(color = "", y = "Rel. No. Cells")

## Joining with `by = join_by(cell_id)`

rel_plt

plot_assemble(
  add_text("(D) Alzheimer plaque density", x = 2.7, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(cowplot::get_legend(spatial_plaq_dens_pl + theme(legend.text = element_text(size = font_size)) + labs(color = "")), x = 43, y = 1, height = 5),
  add_plot(spatial_plaq_dens_pl + guides(color = "none"), x = 0, y = 5, width = 60, height = 61),
  add_plot(umap_plt, x = 52, y = 5, width = 40, height = 61),
  add_text("(E) Expression of \\emph{Jun} depends on the plaque density",
           x = 95, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_text("Pred.~\\emph{Jun} diff.~expr.",
           x = 97.5, y = 6.5, fontsize = 7, vjust = 1, fontface = "plain"),
  add_plot(jun_de_pl, x = 95, y = 4.5, width = 25, height = 61.5),
  add_plot(expr_comparison_pl + guides(color = "none") + theme(axis.text.x = element_blank()), x = 120, y = 5, width = 50, height = 30),
  add_plot(rel_plt + guides(fill = "none"), x = 120, y = 34, width = 50, height = 26),
  add_plot(cowplot::get_legend(expr_comparison_pl), x = 125, y = 60, width = 60, height = 5),
  
  add_graphic("../plots/schematic_elements/alzheimer_circle.pdf", x = -1,  y = -1, units = "mm"),
  
  width = 170, height = 65, units = "mm", show_grid_lines = FALSE,
  latex_support = TRUE, filename = "../plots/continuous_covar_figure_part2.pdf"
)

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## gg[gg1]
## 
## gg[gg2]
## 
## gg[gg3]
## 
## gg[gg4]
## 
## gg[gg5]
## 
## gg[gg6]
## 
## gg[gg7]
## 
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## gg[gg8]
## 
## gg[gg9]
## 
## gg[gg10]

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE